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Abstract 

We present a comprehensive approach to the dynamics of heavy quarks in a quark gluon 
plasma, including the possibility of bound state formation and dissociation. In this ex¬ 
ploratory paper, we restrict ourselves to the case of an Abelian plasma, but the extension 
of the techniques used to the non Abelian case is straightforward. A chain of well defined 
approximations leads eventually to a generalized Langevin equation, where the force and 
the noise terms are determined from a correlation function of the equilibrium plasma, and 
depend explicitly on the configuration of the heavy quarks. We solve the Langevin equation 
for various initial conditions, various numbers of heavy quark-antiquark pairs, and various 
temperatures of the plasma. Results of simulations illustrate various expected phenomena: 
dissociation of bound states as a result of combined effects of screening of the potential 
and collisions with the plasma constituent, formation of bound pairs (recombination) that 
occurs when enough heavy quarks are present in the system. 

Keywords: Heavy Quarks, Quark-Gluon Plasma 


1. Introduction 

Heavy quarks produced in ultra-relativistic heavy ion collisions are interesting for a 
variety of reasons. They are created through hard processes taking place in small space 
time regions, at the very beginning of the collisions, and their abundances remain essentially 
frozen for the entire duration of the collisions. Thus heavy quarks can be used to diagnose 
the properties of the matter they cross before hadronizing. Heavy quarks can make bound 
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states, such as the J/'I' or the T mesons, and the formation of such bound states can 
be strongly affected by the presence of a high temperature quark-gluon plasma. If the 
temperature of such a plasma is high enough, the binding forces can be screened at very 
short distance, hindering bound state formation, as was originally advocated in [1], It was 
also argued that the presence of hot matter in the vicinity of a J/di meson could produce 
an observable mass shift [2], Note that other mechanisms, besides screening, involving in 
particular the collisions with the plasma constituents, can lead to bound state dissociation. 
This is so for instance of the analog of photo-dissociation, namely gluo-dissociation [3, 4, 5] 
(see also [6] for a recent study). 

In situations where heavy quarks are abundantly produced, an excess of bound states 
could occur, due for instance to an enhanced recombination of cc pairs into J/T mesons 
at hadronization. Such a possibility was pointed out early on in Ref. [7] (see also[8]). 
Amusingly, the concern at that time was that the recombination mechanism could spoil 
the proposed signature of quark-gluon plasma, by hiding the expected J/T suppression. 
Recombination was studied systematically using kinetic equations, viewing bound state 
formation and dissociation as a chemical reaction [9, 10]. A more extreme point of view is 
that bound states never truly form in a plasma (we shall come back to this important point 
shortly), but only when matter hadronizes, recombination being then treated as a statistical 
process [11] (see also [12] and [13]). Note that, in contrast to the initial worries, evidence for 
recombination would indicate that heavy quarks roam over long distances through the quark 
gluon plasma before recombining, thereby revealing a rather direct picture of a deconfined 
medium. 

It turns out that the predicted phenomenon of J/T suppression was observed experi¬ 
mentally in the first heavy ion collisions at the SPS, and later at RHIC (the reviews [14, 15] 
include a discussion of experimental results from SPS to RHIC). More recently, evidence 
was obtained at the LHC for sequential dissociation of the T bound states, with the less 
bound states being more suppressed than the most tightly bound ones [16]. Some evidence 
for recombination was also presented by the ALICE experiment [17]. The interpretation of 
the data remains as of today uncertain for a variety of reasons, most of which having to do 
with the production mechanism of the heavy quark bound states (involving issues related to 
structure functions, shadowing, etc), the heavy ion reaction dynamics, etc, all aspects that 
are beyond the scope of the present paper. However, the quality of the recent data, and 
the potential of upcoming experiments, provide strong motivation for further theoretical 
efforts. 

Many investigations concern the fate of the QQ bound state immersed in a quark gluon 
plasma in equilibrium at some temperature T. Such studies where initiated in [18] with 
the determination of the stationary states of a Schrodinger equation, with a temperature 
dependent potential that accounts for screening and the disappearance of the string ten¬ 
sion at high temperature. There has been discussion on the ambiguity in the choice of the 
appropriate potential (free energy versus internal energy) and how to relate it to quan¬ 
tities calculable on the lattice. A review of such potential models can be found in [19]. 
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A somewhat similar line of research concerns the calculation of the spectral functions of 
charmonium states on the lattice. While such calculations have the virtue of being first 
principle calculations within Quantum Chromodynamics (QCD), they suffer from a large 
uncertainty in the reconstruction of the spectral function through the maximum entropy 
technique (for a recent review of lattice calculation at finite temperature, including a dis¬ 
cussion of this issue, see e.g. [20]). Going somewhat beyond the Schrodinger picture, the 
in-medium T-matrix approach [21] can encompass many effects beyond the screening of the 
potential; it can in principle deal with dissociation reactions, changes in thresholds related 
to shifts of masses, coupled channels, etc. 

A recent progress in the direction of a more complete dynamical approach based on first 
principles was initiated in [22]. There one calculates a correlator that is directly related to 
an observable, the rate of dilepton production, and derive the Schrodinger equation obeyed 
by this correlator. A remarkable feature of this equation is that the potential that enters it 
has an imaginary part that reflects the effects of the collisions that the heavy quarks suffer 
with the plasma constituents. The origin of this complex potential was further studied in 
[23], and also in the context of the non relativistic heavy quark effective theory in [24, 25]. 

However, an important issue rarely addressed in the approaches that we have mentioned 
is that the process of the bound state formation is not instantaneous: heavy quarks start 
to interact with the plasma while the correlations that could eventually lead to a bound 
state build up (see e.g. [26] or [27] for early discussions of this issue). This is an important 
feature that should be taken into account when trying to get a complete dynamical picture. 
In short, it is clear that the heavy quarks will suffer collisions with the constituents of the 
surrounding plasma as soon as they are created, and the real issue is whether they will still 
form a bound state when the plasma has cooled down, not whether the bound state will 
“survive”. 

The goal of the present paper is to address this and other issues, by developing a com¬ 
prehensive approach of the entire dynamics of the heavy quarks, including the possibilities 
for bound state formation and dissociation. We shall do that trying to stay as close as 
possible to first principles, and using a chain of well defined (in some cases well controlled) 
approximations. The main objective is to get a global view of the dynamics of heavy 
quarks. As we proceed we shall recover some of the many pieces that have been addressed 
separately in some of the works mentioned above. The approach builds up, extends, and 
to some extent completes, previous works by some of the authors [23, 28]. It is similar in 
spirit to analogous recent efforts using the language of open quantum system [29, 30, 31] 
(see also [32, 33]). We shall be led eventually to formulate, at the end of our approximation 
chain, a generalized Langevin equation, and in that respect, our work bears similarities 
with previous studies using Langevin dynamics for heavy quarks [34, 35, 36]. The present 
work goes beyond such studies by taking into account the dependence of the noise term on 
the configuration of the heavy quarks at each time step. This is an important aspect of the 
dynamics, but it makes the Langevin equation more difficult to solve. In devising suitable 
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techniques to handle it, we were inspired by similar problems in other helds of physics, in 
particular by techniques used in soft matter physics [37]. 

In this exploratory work we focus on the general issue of the formation of bound sates 
of heavy particles in a thermal bath of light particles. The paper is organized as follows. 
In the next section we summarize the general formalism that we are using. We treat the 
heavy particles as non relativistic, and the plasma particles as relativistic. We ignore the 
specihcs of QCD interactions, and for instance the change of their nature (from attractive 
to repulsive) depending on the color state (singlet or octet) the quark-antiquark pairs 
are in. Only Coulomb interactions are retained. The goal of this section is to write the 
probability for a collection of quarks started at some positions at initial time to be found 
at a later time at some other positions. This is formulated in terms of a path integral 
from which we derive an effective theory for the heavy quarks, usually referred to in such 
a context as the “influence functional”. The main approximation in the elimination of the 
plasma degrees of freedom is a weak coupling ansatz that allows us to ignore non linear 
interactions among the gluons. With this approximation, the effective action for the heavy 
quarks takes the simple form of an action quadratic in the charge density of the heavy 
quarks. The whole information about the plasma is contained in a 2-point function. The 
content of the functional is analyzed in section 3 where we study the inhnite mass limit of 
a related object, the correlator of a heavy quark-antiquark pair. In this case, the plasma 
2-point function reduces to a complex potential whose real part describes the screening 
phenomenon, while the imaginary part takes into account the effects of the collisions. The 
inhnite mass limit is the leading order of a systematic approximation, the low frequency 
approximation, that we use to calculate the inhuence functional in the following section. 
The low frequency approximation exploits the fact that the mass of the heavy quark is 
large, and can be viewed as an expansion in terms of the velocity of the heavy particles. 
The leading terms yield a generalized Langevin equation with a multiplicative (position 
dependent) noise. The ingredients of the Langevin equation are related to the real and 
imaginary part of the potential. In Sect. 5 we present results of a set of simulations that 
we have carried out for the case of 2, 10 or 50 pairs of particles. The various aspects of 
heavy quark propagation in a quark gluon plasma are illustrated, including the competing 
phenomena of dissociation and recombination. The last section contains a brief summary. 

2. General setup 

The general problem that we are addressing is that of the dynamics of a collection of 
heavy charged particles in a thermalized bath of light particles with which they interact. 
Although our ultimate goal is to treat QCD interactions, in this exploratory paper we 
focus on the case of Abelian interactions, that is, strictly speaking our picture applies 
directly to electromagnetic interactions. Still we shall use the language of QCD, and call 
the heavy charged particles quarks (with positive unit charge) or antiquarks (with negative 
unit charge). Similarly the plasma in which the heavy particles move will consist of massless 


4 



charged particles, referred to as light quarks, and the photons they exchange will be called 
gluons. We hope this abuse of language will not cause confusion. It is motivated by the 
fact that the simulations that we shall present in this paper involve parameters that are 
adjusted so as to lead to orders of magnitude that are comparable to what we expect for 
heavy quarks in a quark-gluon plasma (in particular we use the strong coupling constant 
as ~ 0.3, not the electromagnetic one a ~ 1/137). Most of the effects that we want to 
study would occur already in Abelian plasmas. Specific effects due to non Abelian color 
interactions will be discussed in forthcoming publications. 

In this section we outline the general formalism that we use. The present approach builds 
on, and extends, that developed by some of us in Refs. [23, 28], and the notation used here 
is close to that of these papers. A related effort was undertaken recently in Refs. [29, 30, 31] 
which directly address the case of non Abelian plasmas. Our goal is to obtain an effective 
theory for the heavy quarks, by eliminating the plasma degrees of freedom. This will be 
achieved by exploiting the fact that the heavy quarks behave as non relativistic particles, 
whose number is conserved as they interact with the quark-gluon plasma, and whose mass 
is large compared to the scales that characterize the plasma dynamics. Ultimately, the 
plasma properties will enter the effective theory only through specific correlation functions 
that are, in leading order, simply related to the potential whose real part describes the 
screening phonemenon, and the imaginary part the effects of collisions. 

The heavy quarks, as we just mentioned, are treated as massive non relativistic parti¬ 
cles. When they thermalize, their typical wavelength, A ~ 1/y/MT with M is the mass and 
T the temperature, is small compared to the inter particle distance of the plasma particles, 
1/T. This suggests that the dynamics of the heavy particles is to a large extent classical, 
and indeed the approximation scheme that we shall present will lead us to a semi-classical 
description. The heavy particles interact among themselves, and with the charged plasma 
constituents via Coulomb interactions. We neglect magnetic interactions, which are sup¬ 
pressed by powers of the velocity, or the inverse mass of the heavy particles. Such restriction 
do not apply a priori to the light particles, but since heavy quarks will not excite magnetic 
modes, we ignore these altogether. Thus the plasma is modeled by massless quarks and 
antiquarks, interacting also with Coulomb interactions. As we shall see, the dynamics of 
the plasma is then characterized by a unique energy scale which is the screening mass rriD- 
In order to get sensible orders of magnitude later in our simulations, we choose this to be 
given by its perturbative value for a two flavor quark-gluon plasma, i.e., = (4/3)5(^T^ 

where g is the gauge coupling. We shall assume throughout this paper that rriD < T M. 

In Coulomb gauge V • A = 0, the Hamiltonian of the system reads 

^ ^ (Pf + Pi) + / (“7^ ^ ^ 

+ \ JJ (2-i) 
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( 2 . 2 ) 


where a* = 707* is a Dirac matrix, and = p + p^ \s the total charge density, with 

N 

p{x) =gY^ [5(* - qj) - 5{x - qj)] , 
i=i 

the charge density of the heavy qnarks and antiqnarks, and 

pp{x) = g (2.3) 

the density of the charged plasma particles (here, charged light qnarks and anitqnarks). 
The plasma is snpposed to be electrically nentral, that is, it contains the same nnmber of 
light qnarks and antiqnarks. It is nsefnl to rewrite the Hamiltonian (2.1) by separating its 
varions contribntions as follows 


H = Hq + Hi+ (2.4) 

with Hq describing the dynamics of the heavy qnarks in the absence of the plasma. 


N 

(Pi+Pi) + ^ // ^xdyp{x)K{x-y)p{y), 
i=i 

Hi the Hamiltonian conpling the heavy qnarks to the plasma charged particles 


Hi = 


jj dxdyp{x)K{x-y)p^{y), 


and Hpi the Hamiltonian of the plasma in the absence of the heavy qnarks 


Hpi = j dx^l^\x) 


(2.5) 


( 2 . 6 ) 


+ m7o^ V’(®) + ^ I I dxdyp^{x)K{x-y)pp{y), (2.7) 


We represent the heavy particles in first qnantization (they are non relativistic particles 
whose nnmber is conserved), while the light particles are represented by the fermion fields 
ip{x) and 'ip^{x). Note that the interaction term in Eq. (2.5) contains contribntions of self 
interactions. Snch terms will not contribnte in the final eqnations that enter onr simnla- 
tions^. We call qj and qj, with j = 1, • • • ,N, the coordinates of, respectively, the heavy 
qnarks and antiqnarks, and pj, pj the corresponding momenta. We denote collectively 
these coordinates by a 2N dimensional vector Q = (qi,--- ,qN,qir ' ^Qn), and often 
refer to Q as a confignration. The last term in Eq. (2.1) is the total Conlomb energy, with 

K{x-y) = --j^ - 1 , -VlK{x -y) = S{x-y). (2.8) 

47 r|a; — y\ 


^They play a role in the real part of the potential to be discussed in Sect. 3. 
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We are interested in the probability P{Q to find the heavy particles in a 
configuration Qj at time tf, given that they are in a configuration at time U. This 
probability may be written in terms of the density matrix of the entire system (the plasma 
and the heavy particles). Let T> be this density matrix. At time ti we assume that the 
heavy quarks have not yet interacted with the plasma, so that the density matrix takes the 
factorized form T>(0) = v'q (8> ^ with P^q = |Qj)(Qj|, a projector on the configuration 

Q,, and is the density matrix of the plasma in thermal equilibrium at 

temperature T = 1//3, with Zpi = exp(/3Tpi) the partition function of the plasma and Fp\ 
its free energy. The density matrix at time t is given by P{t) = e“*'^*T’(0)e*'^*, and the 
looked for probability can be written in the form 

P{Qf,tf\Qi,ti) = Tr {[|Q^)(Qj| (8)1] T>(t)} 

n,m 

^ e-^®™ ... 9 

= |Qi,m)| , (2.9) 

n,m 

where we have set t = tf — ti. In order to trace out the degrees of freedom of the plasma, 
as implied by the formula above, it is convenient to rewrite this expression in terms of path 
integrals. 


We shall do so in steps, in order to identify the main components of the formalism. Let 
us first consider the simple case where the heavy particles interact only with an external 
potential Aq{x), with a Hamiltonian Hi = 9 /o(®)^o(*)- In this case^ 

P{Qf,tf\Qi,ti) = \{Qf,tf\Qi,ti)f . (2.10) 

and the probability amplitude {Qf,tf\Qi,ti) is given by a Feynman path integral 

{Qf,tf\Qi,ti) = (2.11) 

JQi 

where the paths Q{t) satisfy Q{ti) = Qi and Q{tf) = Qf. The actions S'o and 5i are given 
by 


So[Q] 





( 2 . 12 ) 

(2.13) 


^As we shall see shortly this case is relevant for the general discussion. 
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where we have set x = (t, *), and p{x) is the charge density of the heavy particles 


N 

p{x) = '^g[S{x - qj{t) - 5{x - qj{t)] , (2.14) 

i=i 

so that Si can also be written as 


Si[Q,ylo] = -9 



f 

dt 


[MQjit)) - MQjit))] ■ 


(2.15) 


The probability (2.10) can be represented by a very similar formnla, by nsing the closed¬ 
time path formalism. We introdnce a contonr in the complex time plane, as illnstrated in 
Fig. 1, and consider the coordinates {qi,qi} as fnnctions of the complex time tc rnnning 
along the contonr. Alternatively, we may keep time real, bnt dnplicate the coordinates, 
denoting by Qi = {qip,qip} and Q 2 = {qi, 2 ,qi, 2 } the coordinates of the heavy particles 
living respectively on the npper branch (Ci, corresponding to the amplitnde) and the lower 
branch (C 2 , corresponding to the complex conjngate amplitnde) of the contonr. We can 
then write 


P{Qf,tf\Qi,ti) = y (2.16) 

where the actions Sq and Si are given by the formnlae (2.12,2.13) in which the time in¬ 
tegrations are replaced by integrations along the Schwinger-Keldysh contonr. Thns, for 
instance, Sq is given by 


N 

So[Q] = yY^ [ (^1 + 

where the time coordinate P rnns along the contonr C, that is, from ti to tf slightly above 
the real time axis, and retnrns from tf to ti slightly below it. Thns, 

f rtf+iri rU—ir] rtf 

/ dP q] = / dP q] + / dP q] = dt (q|i - q| 2 )> ( 2 - 18 ) 

JC Jti+ir] Jtf—ir} Jt^ 

where in the last step we have dnplicated the coordinate qj{t) as discnssed above. The last 
expression appears natnrally in the action when one mnltiplies the probability amplitnde 
by its complex conjngate in order to bnild the probability (2.16), with qjp{t) labelling the 
path in the amplitnde and qjp{t) the path in the complex conjngate amplitnde. 



Figure 1: The Keldysh contour C, with its different branches. 


It is straightforward to extend the formula (2.16) to include the interactions among the 
heavy particles and with the light plasma constituents. We have (see e.g. [38] or [39]) 

P{Qf,tf\Q„U) = tA) , (2.19) 

where the contour now includes a vertical piece, C 3 corresponding to the thermal average 
of the plasma degrees of freedom at the initial time (i.e., the trace over the equilibrium 
density matrix of the plasma). Accordingly, the fermionic fields in Eq. (2.19) obey anti- 
periodic boundary conditions on C 3 , 'i/j( 0 ,x) = —'0(—i/3, a;), '(/3(0, a;) = —^(—i/3, a;). The 
action 5[Q, V’, V'] is given by 

S[Q,'tjj,'ip] = So[Q] + j d^x ^(a;)( i7^(9^ - m )V^(x) 

p^^^ix)K{x-y)p^^^{y) , ( 2 . 20 ) 

where K{x — y) = 5{tx — ty)K{x — y) represents the (instantaneous) Coulomb interaction, 
and ptot is the total charge density. It is important to stress that the heavy particles do not 
take part in the thermal average, and consequently they do not propagate along the imag- 






inary time sector of the Keldysh contour^. We may take p{t = —ir, x) = 0, with 0 < t < j3. 


The next step consists in eliminating the light fermion field in favor of a Coulomb 
potential Aq. To this end, we use the formal identity^: 


exp 


2 Ptot ■ ■ Ptot 


= AA y DAq exp 


-Aq-K ^ ■ Aq — \Aq- ptot 


( 2 . 21 ) 


where N {det[V2])l is a normalization constant, K ^{x — y) = —5{x — 2 /)Vy, and we 
use a matrix notation to simplify the formulae, e.g.. 


p- K ■ p = JJ^d'^xd^y p{x)K{x,y)p{y). 


( 2 . 22 ) 


When using this identity, it is important to remember that part of the Aq potential is the 
Coulomb field created by the heavy particles. We shall call Aq this contribution, and write 
accordingly Aq = Aq + Aq. By definition, we have 


-V^A^\x) = p{x), A^\x) = j dyK{x - y)p{y). 


(2.23) 


The integration over Aq is then truly an integration over the field Aq, and this field satisfies 
the imaginary time (KMS) periodic boundary condition Ao(0,a;) = Ao(—i/3,®). We could 
take this explicitly into account by performing a shift of integration variables. It is more 
convenient not to do so, provided that we remember that Aq^ plays the role of a constant 
in the functional integral, in particular the result of such integration will depend on Aq . 

Using the identity above, and remembering that ptot = P + P we can perform the 
Gaussian integrals over the light fermion fields 

J exp i Jdx 'ilj{x){ij^dfj_ — m — gj^AQ{x))'tp{x) 

= exp [Tr In — m — py^Ao]] . (2.24) 

The probability (2.19) takes then the form 

P{Qf,tf\Qi,ti) = J J HAo (2.25) 


®Note that we use the notation to denote either a path integrals where the paths are defined on the 
contour, as in DQ, or an ordinary integral, as in dt"^ where the time variable lives on the contour. 
■^We follow closely here the approximation scheme developed in Ref. [28] 
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where Si is given by Eq. (2.13) and 

S 2 [Q, Ao] = -^ j dx{ Ao{x)\7^Ao{x )) - i Tr In [ - m- gj^Aoix) ] . 

The dependence of S 2 on Q arises from the dependence of on the positions of the heavy 

particles. 

It is convenient to rewrite the integral over Aq as the exponential of an effective action, 
the so-called Feynman-Vernon (FV) influence functional [40]: 

gi<I>[Q,Ag'] _ J g-i/cd"‘a;p(x)Ao(x) giS2[Ao]_ (2.26) 

The exponential of the FV functional is the thermal average over the Aq fluctuations of 
the exponential factor that contains the linear interaction pAg between the heavy parti¬ 
cles, with charge density p, and the total Coulomb field Aq. This particular structure is a 
consequence of the fact that the heavy quark is linearly coupled to the total field Aq. So 
far, no approximation has been made (within the present Abelian context). We shall now 
introduce the main approximation of the whole approach, that consists in neglecting the 
non linear self-interactions of the Aq field. 

The action S 2 contains a non-local term describing the coupling between light quarks 
and gluons, as well as the classical Coulomb interaction between the heavy particles. Its 
expansion in powers of Aq gives rise to induced effective couplings to all orders in the 
coupling constant g. In order to be able to compute the influence functional we need to 
introduce some approximations. We do so by retaining only the terms up to quadratic order 
in the coupling g, or equivalently in the field Aq. This is certainly an excellent approximation 
at truly weak coupling, like in electromagnetic plasmas. In the case of QCD, the validity 
of this weak coupling approximation may require further investigation. The main virtue 
of this approximation is to make the path integral over Aq calculable, since it becomes 
Gaussian. The influence functional 4>[(5] becomes 

d'^xd'^y p{x)A^{x - y)p{y) , (2.27) 

where 

^c{x - y) = i{Tc[AQ{x)AQ{y)]) (2.28) 

is the longitudinal gluon propagator (see next section for an explicit expression) defined on 
the contour, and obeying the KMS conditions. Its inverse involves the 1-loop longitudinal 
photon self-energy Hqq, also defined on the contour, 

-A-^{x-y) = K-^{x-y)+lfQQ{x-y), (2.29) 
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where 5^{x — y) = — t^) 5{x — y). 

It is convenient to write the influence functional using the duplicated fields, that is, we 
make the substitution 


^ {Qi{t),Q^{t)) (2.30) 

Ao{t‘^,x) {Ao^i{t,x),AQ^ 2 {t,x)), 

where denotes the curvilinear abscissa parametrizing the Keldysh contour, while t G 
denotes the physical time. The integration over the physical time is always from ti 
to tf . With this notation A becomes a matrix, 

^ab{tx - ty) = - t^) with t^eCa,t^eCb, a, 6 = 1,2. (2.31) 

We have, explicitly, 

Aii(a:, 2 /) = i(r^ [Ao,i(a:)vlo,i(y) ]) = i{T [Ao{x)Ao{y)]) = A{x,y), 

A2i{x,y) = i{T^ [Ao,2(a;)^o,i(y)]) = K^o{x)Ao{y)) = iA>(x,y), 

Ai 2 {x,y) = i{T^ [2lo,i(a;)ylo,2(y)]) = i(^o(y)^o(a:)) = iA<(x,y), 

A 22 {x -y) = i{ [^o, 2 (a:^) 2 lo, 2 (y) ]) = i( ^ [^o{x)Ao{y) ]) = A(x, y). (2.32) 

where T and T denote respectively the time ordering and anti ordering, and (• • •) is the 
thermal average. Using this notation, the phase in Eq. (2.27) becomes 

^[Q] = \ [ / dty f dxdy (-1)“+Va(4, ®) Aab{tx - ty, x - y) Pb{ty,y). 

J J t'i J 

(2.33) 

Note that we integrate all times the same way, i.e., on the real time axis from U to tf, so 
that the off diagonal terms pick up a minus sign. Note also that in the right hand sides of 
Eqs. (2.32), we have introduced the notation A (without subscripts) to denote the time- 
ordered real time propagator. Other useful relations are A = A^ -|- iA^, A = — A"^ -|- iA^, 
where A^ and A"^ denote respectively the retarded and the advanced propagators. 

At this point the probability (2.16) is written as the following path integral 

P{Qf,tf\Qi,ti) = I^DQ (2.34) 

with d>[Q] given by Eq. (2.33) above. The plasma degrees of freedom have been eliminated, 
the plasma properties entering the calculation of [Q] solely through the contour propagator 
Aab{tx ~ ty,x — y) that plays the role of an effective interaction between the heavy quarks. 
The problem of calculating the probability P{Qf,tf\Q^,ti) has been reduced to that of 
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calculating an “ordinary” Feynman path integral. This remains however a difficult task, 
in particular when many heavy QQ pairs are present in the system, and we shall shortly 
proceed with further approximations. Before we do that, we shall consider in the next 
section a situation where the influence functional can be calculated exactly: this is the case 
where a single, infinitely massive, QQ pair is present in the system. 


3. The influence functional and the complex potential 

We consider in this subsection the case of a single heavy QQ pair. This will allow us 
to relate the influence functional to the complex potential that was first identified in this 
context in Ref. [22]. We denote here the coordinates of these heavy particles by Q = {r, r}, 
and we consider the correlator: 


G^{tf,Qf\ti,Qi) = {i^Q{tf,rf)i;Q{tf,rf)ipl^{ti,ri)ip^^{ti,ri)), (3.35) 

where the angular brackets represent the thermal average over the plasma particles (being 
understood that the states of the plasma do not contain any heavy quarks). This object 
enters directly the calculation of the heavy quarkonium spectral function, and for instance 
the calculation of dilepton emission rate [41]. As was shown in [23], under the same ap¬ 
proximations as those done presently, this correlator is given by 

G>itf,Qf\ti,Qi)= (3.36) 

JQ, 

where Q lives on the upper part of the contour. The influence functional has the same 
form as in Eq. (2.33), that is. 


<h[Q] 



dx dy p{tx, x) A{tx -ty,x- y) p{ty,y) , 


(3.37) 


but now all times are on the upper part of the contour, and here A = An is the real 
time, time-ordered, propagator (see Eqs. (2.32)). Recall that the density is p{tx,x) = 
g \5{x — r{tx)) — S{x — r(ta,))], so that the influence functional can be written as ^qq + 

+ ^QQ’ 


<^00 [Q] 




<^00 [Q] 



u 




dt^ A(tx 
dt^ A(tx 


/ dty A{tx 


tyj 'f'itx) 

tyi 'f'{tx) 


r{ty)), 

rity)), 

-r{ty)). 


(3.38) 

(3.39) 

(3.40) 
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In the last line we have used the fact that A{t, x) is in fact a function of \t\ and \x\ in order 
to rewrite A{tx — ty,f{tx) — ^{ty)) as A{ty — tx,r{ty) — f{tx)), which coincides with the 
term already written after the interchange of the integration variables tx and ty. 


The calculation is particularly simple in the infinite mass limit, where the paths are 
trivial (since infinitely heavy quarks do not move). In this case, the influence functional 
takes the form 

^[Q] = -9^ f dtx f dty [A(4 -ty,r -r) - A{tx - ty, 0)] . (3.41) 

t/ J 

At this point, it is convenient to express A{tx —ty,r — f) in terms of its Fourier transform 
A{tx-ty,r -r) = y - r). (3.42) 


This allows us to perform the time integrations 



dtye 


LO^ 


(1 - cos{uj{tf - ti)) , 


and obtain, after a further Fourier transform of the coordinates. 


(3.43) 


m ] = 25' 


do; 


dfc 1 

(27r)^ 


cos{ujt) 




A{co,k) , 


(3.44) 


with t = tf — ti- 

We are interested in the evolution of the heavy quarks on time scales that are large 
compared to the time scale that characterizes the dynamics of the plasma, and which is 
controlled by the inverse of the Debye mass, mo- It is then useful to consider the large 
time limit of the expression above. This is easily obtained with the help of the relation 
(1 — cos(cjt))/a;^ ~ TTt6{Lj) valid as t —)■ 00 (i.e., t 1/w). We get 


l-IQI == sHtf - *•) / (^ (1 - e'""-"’) A(0. fc). 


(3.45) 


Thus, at large time, the influence functional is dominated by the zero frequency part of the 
response function of the plasma. As an alternative to the calculation done above, we could 
change integration variables, tx,ty —)• {tx + ty)j2,tx — ty, and observe that when tf — ti is 
large (compared to on can integrate freely over tx — ty, which filters out the zero 

frequency component of the response. 

By considering the equation of motion for the correlator (3.35) at large time, and its cor¬ 
responding expression in terms of the influence functional, we interpret, following previous 
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works, the coefficient oi t = tf —ti in the influence functional as a complex potential. That 
is, we set e‘^ = More precisely, remembering that A(0,r) = A^(0,r) + iA^(0,r), 

we set 

/ flk 

— e^'^-^A^{u: = 0,k), (3.46) 

/ (ik 

^e‘^-A<(0,fc), (3.47) 

so that 


V(r) = -g‘^[V{r) - Ken(O)] -ig\W{r) - fT(0)]. (3.48) 


where the minus sign in front of g^ is the same as in Eq. (3.40) and reflects the fact that 
the r dependence of the potential describes interation between heavy quarks with opposite 
charges. At this point, we identify A with the real-time gluon propagator in Fourier space, 
at zero frequency. In the hard thermal loop approximation [42], a suitable approximation 
in the present context, this is given by (see e.g. [23]) 


Dl{u: = 0,k) 


— 1 . 

k'^ + m‘j^~^ |fc|(fe^-|-m|,)2 ’ 


(3.49) 


which allows us to get an explicit expression for V(r), a function of r = \r — r\. The 
calculation of the integrals in Eqs. (3.46) and (3.47) yields 


V(r) 


dvr 


rriD 


g^e-mor g 2 T 


(3.50) 


where the first term is a self energy contribution, I4en(0), that has been estimated by sub¬ 
tracting the corresponding vacuum part, thereby leaving the following thermal contribution 


/ (- 4 ) = -T- ( 3 - 51 ) 

Jq\q^ + m% q^J dvr 

Note that, as expected, the real part of the potential between the quark and the anti-quark 
is attractive and screened. The imaginary part of the potential originates from the collisions 
between the light fermions of the hot medium and the heavy quarks. In fact, the quantity 


g^T r Sk nmjjT 

J (27r)3|fc|(fc2 + ^2^)2 


(3.52) 


is the rate of collisions between one heavy quark and the light quarks of the plasma. It may 
be identified with the damping factor associated with the propagation of one heavy quark 
in the plasma. The function 


roo 

4 >{x) = 2 / dz 

Jo 


(z2 + l)2 _ 


1 - 


sm Urc 


zx 


(3.53) 
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which vanishes for x = 0 and increases monotononsly, approaching nnity as x—7-+00. Thns, 
the collisional damping rate is most important when the heavy qnarks are far apart: when 
this is so, the damping factor associated with the propagation of the heavy qnark pair is jnst 
twice the damping factor of a single heavy qnark. When the heavy qnark gets closer to the 
heavy qnark, interference occnrs that gradnally snppresses the effect of collisions. When the 
QQ separation vanishes, this interference is completely destrnctive and kills the imaginary 
part: this is becanse when the QQ separation is too small, the charge of the individnal 
heavy qnarks are not resolved by the typical flnctnation of the electric potential of the 
plasma. The QQ pair behaves then as a small, nentral, electric dipole, which propagates in 
the plasma withont interacting. Note that the behavior of the fnnction 4>{x) at small x is 
singnlar: 4>{x) is continnons as x —)■ 0, bnt it does not have a simple Taylor expansion, as 
can be seen from the logarithmic divergence of the coefficient of the term of order x^. We 
come back to this issne in Sect. (5.1). 

Before we leave this section we shonld emphasize an important difference between the 
calcnlation that we have jnst presented of the correlator (3.36), and that of the probability 
(2.34). Snperficially, these qnantities differ solely by the contonr involved in the integration 
of the inflnence fnnctional. In fact this change of contonr is not innocent, and the two 
qnantities are deeply different. The correlator (3.35) may be interpreted as a probability 
amplitnde to find the QQ pairs in confignration Qj at time tf together with the plasma in 
the same state as it is at time ti, irrespective of what that state is. Becanse, dnring their 
propagations, the heavy qnarks mix with complicated confignrations that involve plasma 
constitnents, the amplitnde decays with increasing time, and this even when the qnarks 
are infinitesimally heavy. This is the origin of the imaginary part of the potential, an 
imaginary part that wonld affect also the analog of the correlator (3.35) for a single qnark 
[28]. The probability (2.34) addresses another qnestion, namely the probability to find the 
heavy qnarks in the confignration Qj, irrespective of the state of the plasma. In the limit 
of infinitely heavy qnarks, one expects this probability to be proportional to S{Qj — Qj) 
and this is indeed what we shall verify in the next section. 

4. Low frequency approximation and generalized Langevin equation 

We now retnrn to the calcnlation of the probability (2.34), and introdnce an approx¬ 
imation, the low freqnency approximation, that allows ns to go beyond the infinite mass 
limit that we have considered in the previons snbsection. Still, as we shall see, in this ap¬ 
proximation, only the basic qnantities that appear in the infinite mass limit will be needed, 
namely the real an the imaginary part of the potential. 

4-1. The low frequency approximation 

The approximation relies on the fact that the dynamics of the heavy fermions is mnch 
slower than the dynamics of the light fermions of the medinm. As we have recalled, the 
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typical frequency in, for instance, the time-ordered propagator^ is mo- Now, during 

a time t, the heavy quark moves a typical distance y/tjM. For t ~ this is a small 

distance compared to the size of the screening cloud, ~ mriD ~ \/itid/M <C 

1 . Thus, over a time scale characteristic of the plasma collective dynamics, the heavy 
quark positions are almost frozen (they are completely frozen in the limit M —)• oo). Said 
differently, the plasma dynamics looks very fast to the heavy quarks, and their interactions 
with plasma constituents are essentially instantaneous. In order to observe a substantial 
motion of the heavy quarks, we need to look at the system over time scales that are large 
compared to In the calculation of the influence functional, we need therefore to allow 

for tx — ty ^ equivalently, in Fourier space, typical frequencies uj mo- To get 

a systematic expansion, one expands A(a;) in powers of oj around cj = 0. In leading order 
this yields 

A(4 -ty) = J [A(a; = 0) + i^A'{co = 0)] 

- S{tx - ty)A{u) = 0)+ i-^5{tx - ty)A'{uj = 0). (4.54) 

dtx 

This expansion shows indeed that the heavy quarks interact with the medium via effectively 
instantaneous interactions. One recognizes in the first term of the expansion the infinite 
mass limit. The corrections implied by the second term will involve the velocities of the 
heavy quarks, as we shall see shortly. 

We have identified in the previous subsection the zero frequency part of the time or¬ 
dered propagator to the complex potential. Because of the relation obeyed by the various 
components of the propagator (see for instance [43]), the derivative term A'(a; = 0) is 
simply related to the imaginary part of the potential, as we now show. Indeed, the time or¬ 
dered propagator can be written as A(ti;) = A^(w) -|-iA^(w), where A^(w) is the retarded 
propagator and A^(a;) has been defined above. The latter is related to the other function 
A^{lo) by the KMS relation, A>(a;) = e^‘^A^ (cj), and the two functions allow us to re¬ 
construct the spectral density p{uj) = A^(a;) — A^(a;). From the last two equations, one 
easily establishes that A^(a;) = N{uj)p{uj), with N{uj) = l/(e^‘^ — 1). From this relation, 
and using the fact that the spectral function is an odd function of w, it is easy to show that 
A^(—cj) = A^(a;), so that, in particular, A^(0) = A^(0). It follows then easily that 


dA> 

> 

A 

dA< 

dtc 

a;=0 dw 

<..=0 ’ dw 



(4.55) 


Furthermore, it is easily shown using the spectral representation of the retarded function, 
and again the fact that the spectral density is an odd function of io, that dA'^(a;)/da;|^_g = 


®Since the spatial coordinates, or the three momenta, play no role in this discussion, we omit them to 
simplify the notation and denote the propagator A(cj, r) or A(cu, k) simply by A{uj). 
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0. Therefore 


dA(a;) . dA<(a;) 
dw ^^0 da; 


a ;=0 




(4.56) 


Note finally that A 22 (w) = —A^{uj) + iA^(a;), where A^(u)) denotes the advanced prop¬ 
agator. At zero frequency, A^{aj = 0) = A^{aj = 0). Thanks to the relation (4.56), the 
expression (4.54) involves no new ingredient beyond the real and the imaginary part of the 
potential introduced in the previous section. 


After this preparation, we can now calculate the influence functional in the low frequency 
approximation. To do so, we use Eq. (4.54), the relations Aii(O) = A^(0) -|- iA^(O), 
Ai 2(0) = iA<(0), A2 i( 0) = iA>(0) = iA<(0), and A22(0) = -A^(0) +iA<(0 = -A^(0) + 
iA^(O), together with the definitions (3.46) and (3.47). We write the influence functional 
as <h[Q] = ^qq[Q] -|- ^ straightforward calculation then yields 

[Q] = y Z] / - qi,2) - V (Qj -1 - qi,i) 

i,j = l 

-iW{qj^2 - qi,2) - iW^(qi,i - qi,i) + 2[W{qj^i - qi^2) 

+ ^{qi,2 + qj,i)--^W{qj^i-qi^2) ■ (4.57) 

2 Clqi,2 

and similarly for ^qq [Q] with the substitution {qj} —)■ {q*}. The mixed fermion-antifermion 
contribution reads 


/•*/ r 

$ _[Q] = -5^ Z / dt v{qj^2 - qi,2) -y{qj,i - qi,i 

y y / +. L 


i,j=l •' 


-iVE(qj,2 - qi,2) - '^W{qj,i - q^q) -F iVF(qj,i - qj^2) + iW^(Qi,2 - qi,i) 


(4.58) 


Note that in the infinite mass limit, we can identify the coordinates of the heavy quarks 
in the upper branch of the contour with the corresponding ones in the lower branch, e.g., 
qjp = qj^ 2 - Furthermore, in this limit, we can ignore the velocity q. it is then easily verified 
that in this situation the influence functional vanishes identically. And this is as it should. 
We have then P{Qptf\Q^,ti) = S{Qj — Qi). (See also the discussion at the end of the 
previous section.) 


4-2. Generalized Langevin equation 

We shall now use the expressions that we have obtained for the influence functional 
<1> in the low frequency approximation in order to perform a further approximation that 
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will lead us to a reformulation in terms of a generalized Langevin equation. This last 
approximation exploits the fact that the trajectories of a heavy particle in the amplitude 
do not differ much from that in the complex conjugate amplitude. This suggests to perform 
the following change of variables 

= yi = (iii-qi2, (4-59) 

(and similarly for the antiquarks coordinates), and to expand the influence functional in 
powers of the small deviations and y^. In order to motivate this expansion, we note that, 
after an integration by parts, the exponential of the free action takes the form 


exp 


^ rtf 

-i ^ X] / ’ 

i=l 


(4.60) 


and similarly for the antiquarks. The dominant contribution to the path integral comes from 
the region where the phase in Eq. (4.60) is small or at most of order unity. We can estimate 
the integral as dt(fj • ~ y^TjM\y^\, where y^T/M is the thermal velocity of the 

particle. The condition that the phase be small is then that lyj be small, |yj I < 

We then proceed to the expansion of the influence functional to second order in y^. 
The details of this expansion are given in Appendix Appendix A, and we report here the 
result. We collect the coordinates {ri,fi,y^,y^} into 2N dimensional vectors®, as we did 
earlier for Q, e.g. R = (ri, • • • r^r, fi • • • f^r). As a result of the expansion, we can write the 
probability ^(Rjj | Rj ,ti) as follows 


where 


rR,/ /*Y/=0 

DR I DY exp 




'Yi=0 


ftf 


dtC{R,Y) 


(4.61) 


£(R,Y) 


^-i Y • (^MR + M7(R) • R - F(R)) 


1 

2 


Y A(R) Y 


(4.62) 


We have Yj = 0 = Yj because the coordinates and qi ^2 of the heavy particles coincide 
at the ends of the Schwinger-Keldysh contour. 

The 2A^-dimensional vector F(R) represents the forces between the heavy particles. It 
is given in terms of the gradient of the potential V (r) as follows 

TV / VViri - r,) - VViri - r,) \ 

F,,(R) = VE (4.63) 

j=i \ VV{fi - fj) - VV{fi - Tj) j 


®In fact these vectors have 2N x 3 components since for instance ri is a three component vector. We 
do not indicate explicitly these components in order to alleviate the notation. Similarly, for each pair of 
vectors labelled by i and j, say ri and rj, Hapivi — rj) is a 3 x 3 matrix mixing the components of the 
corresponding vectors. 
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where i = 1,..., A^, and the primed index i' runs from 1 to 2N, with i = i' for i' < N { 
first lin e of (4.63)), i = i' — N for i' > N (second line of (4.63)). The first line of Eq. (4.63) 
represents the force exerted by all the heavy quarks and antiquarks on the heavy quark 
at position r*, whereas the second line is the corresponding force exerted on the heavy 
antiquark at position fj. 

The {2N X 2A^)-dimensional matrix 7(R) represents the friction exerted by the medium 
on the heavy particles. Its expression involves the Hessian matrix Ti of the function W, the 
imaginary part of the potential, and reads 


Ti'i' (^) — 


9 


2MT 


nivi - Tj) -nivi - Tj) 
-n{ri - Tj) niri - Tj) 


T-La/sir) 


(9IT(r) 
dradrp ’ 


(4.64) 


where the primed indices i\j' = 1,..., 2N are related to the unprimed ones, respectively 
i and j, as indicated above. The Greek indices a,l3,'y label the cartesian coordinates of 
r. The matrix 7 is symmetric and real (hence diagonalizable with real eigenvalues^). This 
follows from the fact that, for instance, 'H{ri — fj) = ^(fj — r^), and the fact that the 3x3 
matrix ^Q,^(r), being a Hessian matrix, is symmetric. 

Finally, the matrices 7 and A in Eq. (4.62) obey Einstein’s relation 


A(R) = 2MT7(R). 


(4.65) 


In the Appendix Appendix B we show that the probability (4.61) can be generated by 
the following generalized Langevin equation [44] 


M R = -M 7 (R) • R + F(R) + ^(R, t ), 


(4.66) 


with a space dependent (also referred to as multiplicative) white noise ^(R, t): 


(^j/(R, t) ) = 0, (^fc'(R, t) <^(i - t') . 


(4.67) 


The fact that the friction (and hence the noise) depends explicitly on the configuration of 
the heavy quarks is what makes this Langevin equation distinct from what has been done so 
far in the context of heavy quark dynamics. The mathematical subtleties of such Langevin 
equations with multiplicative noise are recalled in the Appendix Appendix B. Let us just 
mention here that the present equation, with its explicit inertia term, does not suffer from 
discretization ambiguities, and we have used the Ito prescription to solve it numerically (see 
Appendices Appendix B and Appendix C for details). 

In order to get a first orientation as to the effect of the spatial dependence of the noise, 
we consider in the next subsection the simple case of a single pair of heavy particles, one 


^We shall see that the eigenvalues are also strictly positive, which is physically expected for a matrix 
representing a friction term. 
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heavy quark and one heavy antiquark, for which analytical results can easily be obtained. 
This will also be used to introduce the notion of a bound state in this classical setting, 
and how such a bound state evolves when it is in contact with a thermal bath at various 
temperatures. 

4-3. Langevin equation for a single heavy quark antiquark pair 

When a single heavy quark antiquark pair is present in the system, the generalized 
Langevin equation (4.66) takes the form 

B ( 7 ^ 

M f + (^(0) r — 'H{s)r) — g"^ W (s) = t) 

M f ^ {n{ 0 ) f - n{s) r) + g^ VV{s) = f{s, t) (4.68) 

where we have set s = r — f, with r and r denoting the coordinates of the quark and the 

antiquark, respectively, and the correlators of the noise are given by® 

( Ca(s, t) C/sis, t')) = { fais, t) fg{s, t') ) = g'^ n{0) 5ag5{t - t') 

( ?a(s, t) t') ) = -g"^ Tiagis) S{t - t') . (4.69) 

The two Langevin equations are correlated through the force terms, as well as the fric¬ 
tion terms which depend explicitly of the distance between the two heavy particles. It is 
convenient to write these equations in terms of relative (s = r — f) and center of mass 
(p = (r-|-r)/2) coordinates. By taking the sum and differences of the two equations above, 
we get 

Mp+^ |W{0) - n{e)] p = + 

Ms + ^ [W(0) + «(«)]-2/vr(3)=j(s,t)-«(<i,t). (4.70) 

Note that only s is sensitive to the (attractive) force between the quark and the antiquark. 
The center of mass of the pair just follows a random walk and is subjected to a drag 
force and a random force. When the size of the pair exceeds the Debye radius, i.e., when 
svriD ;§> 1, T-L{s) ~ 0, and the noises ^ and ^ become uncorrelated. Using the fact that 

g^n{0)^g = 2MTj5^g (4.71) 

is diagonal, with 7 constant, we can then rewrite the equation for p as a standard Langevin 
equation for a particle of mass 2M and drag force 7 . In fact, in the limit smo 1 the 

^Recall that 4^ is a 3 x 3 matrix, and that ^ and f are here three dimensional vectors. 
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two equations (4.68) decouple. At large time, |t — to| the mean square displacement 

of r{t) (and similarly for f(t)) follows then the law of diffusion, 

{{r{t)-r{to)f) = 6V\t-to\, V = (4.72) 

where T) is the diffusion coefficient. 

In the opposite situation where sttid 1, the friction term cancels in the equation 
for p: this is because in this case the quark and the antiquark form an electric dipole of 
very small size that propagates in the plasma as a color neutral particle of mass 2M, and 
hence does not interact with the plasma (one can easily verify that the contributions of the 
random forces also cancel, as they should). In the same limit of small size we can expand 
the potential, V^(r) ~ 1/(0) + (l/2)A:r^, with k = d?V /and rewrite the equation 
for the relative motion as 

+ = C(t), (4.73) 

with ( ((j(f) ) = 0 and 

( at) Ci(t') ) = - t'). (4.74) 


The heavy quark pair behaves then as a harmonic oscillator coupled to a thermal bath. This 
is as close as we can get to the notion of a bound state in this classical picture. Assuming 
that the expansion of the potential to quadratic order remains valid at large time, the mean 
square displacement of s{t) will eventually reach its value in thermal equilibrium, given by 
the equipartition theorem: 



(4.75) 


This formula indicates that the radius of the QQ pair increases with the temperature^*^. 
As the temperature increases, the radius becomes eventually too large for the harmonic 
approximation to the potential to remain meaningful. In fact when this happens, the 
potential becomes essentially flat, indicating that no force maintains the QQ pair together: 
the bound state dissociates. These qualitative predictions will be made more quantitative 
in the next section. 


®We discuss in the next section how to regulate the Coulomb potential so as to give meaning to this 
expansion. 

^®We shall see in the next section that the explicit linear dependence on T of the numerator is in fact 
amplified by the decrease of the coupling constant, and to a less extent that of the spring constant k, with 
increasing temperature 
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5. Numerical results 


We now present results of simulations of systems containing a given number N of heavy 
quark antiquark pairs in a quark-gluon plasma at temperature T. We discuss first the case 
of a single pair, N = 1, and follow its fate for various temperatures, thereby turning the 
considerations of the previous subsection into a more quantitative discussion. Then we 
turn to the case of many pairs (up to N = 50), where, in addition to the phenomenon of 
dissociation that occurs for a single pair, the formation of new bound states through the 
process of recombination is also possible. 


The parameters in the problem are the mass M of the heavy quarks, the temperature 
T of the plasma, and the gauge coupling g. In thermodynamical calculations, the latter 
quantity depends on the temperature and is commonly chosen to be the running coupling 
at a scale ~ 2'itT. Although this is not crucial in the present work, we take into account 
this running of the coupling with the following simple relation taken from Ref. [45] 


a, = 


dvr 


Ols{Tc 


1 + C In ( ^ 


C = 0.760, Tc = 160 MeV, as{Tc) = 0.5. 


(5.76) 


The Debye mass is approximated by its perturbative expression for a two flavor quark gluon 
plasma, | . With the running coupling given above we have uid 460 MeV 

for T = Tc- The coupling of the heavy quark to the plasma constituents involve an extra 
color factor Cp = 4/3 which is ignored. Finally we shall consider charm and bottom heavy 
quarks, whose masses are taken to be respectively Me = 1.4 GeV and Mf, = 4.2 GeV. 
Again, we emphasize that all these numbers, as well as all those which follow in this entire 
section, are meant to provide reasonable orders of magnitude, in line with those expected 
for quarkonia in a quark-gluon plasma; but we are not attempting to develop here a precise 
phenomenology. 

In addition to the physical parameters that we have just discussed, we need to specify 
another one, a cutoff A, whose role is to control the short distance behavior of the real and 
imaginary parts of the heavy quark potential. This requires more discussion, and is the 
object of the next subsection. 


5.1. Estimation of the cut-off 

Before we can use the Langevin equation derived in the previous section, we need indeed 
to cure two problems associated with the short distance behavior of the complex potential. 


Gonsider first the real part V (r) . This is given by the screened Goulomb potential, which 
behaves as 1/r at short distance. This poses a well known problem in classical simulations. 
One way to see it is to notice that the classical distribution, that the Langevin equation 
eventually leads to, ~ is singular at small r for the attractive Goulomb potential. 


23 




X (fill) 


Figure 2: The potential energy of a QQ pair as a function of the Q-Q distance x for three different 
temperatures, and calculated with a cutoff A = 4 GeV. Much of the temperature dependence of the potential 
at small distance {x ~ 0.05/m) can be attributed to that of the running coupling. The temperature 
dependence of the screening mass mo ~ T affects the potential in the intermediate range {x ~ 0.4/m). 
The oscillations at intermediate and large distances are an artifact of the finite cutoff. 


This would lead to an infinite probability for two particles to be close together. We may 
also observe that when two particles come to close to each other, their relative kinetic 
energy becomes big, and this violates the conditions of validity of the approximations used 
in Sect. 4 when deriving the classical equations. Note that this is a problem that arises 
only in the classical treatment of the Coulomb interaction through the Langevin equation; 
it would not occur if we were to solve the corresponding Schrodinger equation. A simple 
way out is to add a repulsive “quantum correction” in the form h?/2Mr'^, as originally 
proposed by Kelbg [46]. Many refinements of this procedure have been studied (in the 
present context, see for instance [47] and references therein). In this exploratory work, 
we find it sufficient to turn off the force at short distance, as was done for instance in 
[35]. We do so here by introducing a finite cutoff in the integral Eq. (3.46) that yields the 
screened Coulomb potential. The resulting potential is displayed in Fig. 2. Note that when 
calculated with this prescription the value of the potential at the origin, 1^(0), depends 
linearly on the cutoff, g'^V{0) ~ {2as/Tr)A. Therefore the cutoff A cannot be chosen too 
small otherwise the potential will not be deep enough to sustain bound states of the bottom 
quarks. Taking this into account, as well as further consideration to be presented shortly, 
we have settled for a value A = 4GeV, and this is the value with which the plots in Fig. 2 
have been done. The temperature dependence that is seen in Fig. 2 arises mainly from the 
temperature dependence of the coupling constant, according to Eq. (5.76). 

The presence of the cutoff makes the potential regular at short distance. One can then 
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expand it around the origin and find the spring constant k introduced in Sect. 4.3. We get 


k I ( K A \ 

^ -^ arctan- . (5.77) 

rri'jj DTT^ yimjy mo mo J 

As mentioned in Sect. 4.3, the bound state will dissociate when the size, as measured by 
(s^) = 3T/ {g^k) becomes of the order of the Debye radius, m~j^. Defining the corresponding 
temperature as Td, we get (when A » m^) 


Td 




k 


(5.78) 


For A = 4 GeV, rriD = 0.5 Gev, this yields Tjy = 320 MeV, a reasonable order of magnitude. 
This provides another argument in favor of a not too small cutoff. 



A (GeV) 


Figure 3: The cutoff-dependence of the diffusion coefficient V — T/My (see Eq. (4.72)), multiplied by 2tiT, 
for different values of the temperature. Note that the differences between the three curves corresponding 
to different temperatures is largely due to the variation of the coupling constant, according to Eq. (5.76). 


The second reason why we need a cutoff is that the second derivative of the imaginary 
part of the potential, that enters the definition of the friction, is divergent. This was already 
mentioned at the end of Sect. 3. The problem here is of a different nature as that of the real 
part. It reflects the fact that the hard thermal loop approximation used in the calculation 
of the imaginary part of the potential involves kinematical approximations that cease to 
be valid whenever large momentum exchanges are involved. Again the divergence can be 
controlled by a cutoff, which, here, would be naturally of the order of the temperature. In 
fact, we shall proceed as for the real part of the potential, and simply limit the momentum 
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integral in Eq. (3.47) to values lower than A. Note that the values of A that are needed 
for V and W are a priori unrelated to each other. However, for simplicity and in order 
to avoid the proliferation of irrelevant parameters, we have performed calculations with a 
common value for A, independent of the temperature. It turns out that the drag coefficient 
and the diffusion constant depend only mildly on A around the value A = 4 GeV that we 
have adopted (see Fig. 3). 

From the second derivative of W we can calculate the drag coefficient, according to 
Eq. (4.71) and we get 


_ mfy / / A^ \ A^ \ 

24 TT M \ \ ) A^ + ) ’ 


(5.79) 


To within a color factor Cp that we ignored, and with the specific choice A = T, this 
expression agrees with that obtained in Ref. [34] in the leading logarithm approximation. 
The diffusion constant T> = T/{M 7 ) is plotted in Fig. 3 as the dimensionless combination 
{2ttT)V-. 


V ■ 2ttT = 


9 

4a2 




A2 


+ 1 


-1 


(5.80) 


One sees that in the region A ~ 4 GeV, the diffusion constant depends indeed weakly on 
the value of A. Furthermore, for this value, 2'kT'D ss 2.7 for T = 160 MeV, or 7 ss 0.2 fm“^. 
These values are of the order of magnitudes of those used in phenomenological studies [34] 
(see also [48] for more recent estimates). 


Now that we have adjusted all the parameters, we can start exploring the main features 
of the dynamics of the heavy quarks in a plasma, as predicted by the generalized Langevin 
equation (4.66). The details of the numerical method that we use to solve this equation are 
given in Appendix Appendix G. 

5.2. One heavy quark-antiquark pair 

Our first set of results concerns the evolution of a heavy QQ pair immersed in a uniform 
quark-gluon plasma in thermal equilibrium at temperature T. The pair is prepared so that 
it corresponds initially to a bound state with a given size and binding energy. One first 
generates a sample of pairs, with the following procedure: The distance between the quark 
and the antiquark is chosen randomly between 0 and the Debye radius rp = ni^. The 
relative initial velocity of the quark and the antiquark is taken from a Maxwell distribution 
centered at the average value of typical quarkonia relative velocities (e.g. Vq ~ 0.3 for 
charmonium [49]). Then we select from this sample the pairs that can be associated with 
specific bound states according to criteria that will be specified shortly. We then simulate 
the evolution of the pair using the Langevin equation (4.66) that was derived in Sect. 4.3. 
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As a first check of the Langevin dynamics, we consider a cc pair at a temperature 
T = 200 MeV. At this temperature, and for the parameters that we have chosen, all cc 
bound states eventually dissociate in the plasma. This is what the plot on the left of Fig. 4 
indeed shows. After an initial transient period of time, the two consituents of the pair 
follow independent Brownian motions, with the average distance squared growing linearly 
with time, in agreement with the analytical result, Eq. (4.72). Moreover, the right panel 
of Fig. 4 shows that the constituents indeed thermalize, the energy per quark reaching the 
value (3/2)T, in agreement with the equipartition theorem. 



0 10 20 30 40 50 

t (fm/c) 



_I_I_I_I_ 

0 10 20 30 40 50 

t (fm/c) 


Figure 4: On the left: Average c-c distance squared as a function of time. This follows the predicted long¬ 
time Brownian behavior with the diffusion constant given by T> ■ 2'kT ~ 4. On the right: Average energy 
(in 200 MeV units) per quark as a function of time compared with the energy at equilibrium (horizontal 
line). In both graphs we used T = 200 MeV and A = 4 GeV. Statistical errors are to small to be plotted. 


However, the very long time, where the heavy quarks eventually thermalize with the 
surrounding plasma, is not our main concern here. We want to understand the dynamics 
over shorter time scales, in particular because the plasma produced in a nucleus-nucleus 
collisions has a hnite lifetime. To be specihc, we shall take this lifetime to be Tqgp ~ 10 fm/c, 
and accordingly our main focus will be to understand the dynamics of the heavy quarks 
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over such a typical time scale. We shall also differentiate between different charmonium 
states, J/'h (IS), Xc (IP) and 'h' (2S), but consider a single bottomonium state which 
we shall refer to as the T. A word of explanation is needed here regarding what we mean 
by bound states. Within the classical simulation using the Langevin equation, this refers 
to the following procedure. At the beginning of the simulation we calculate the binding 
energy of a pair in its center of mass frame^^ and select the pairs according to the values of 
their initial radius tq, and their binding energy AE. Depending on these values, we call a 
pair by the name of the closest bound state it would correspond to in a complete quantum 
treatment. The specific criteria that we use to attribute a charmonium state to a given pair 
are the following 

• T' : 0 < AE' < 100 MeV and tq > 0.35 fm, 

• Xc ■ 100 < AE < 300 MeV and vq > 0.25 fm, 

• J/T : AE > 550 MeV and ro > 0.10 fm, 

• T : AE > 700 MeV. 

For the bottomonium, as already mentioned, we do not attempt to discriminate between 
the various bound states, and the requirement of a large binding energy automatically 
selects small sizes. In the case of charmonia, the constraint on the radius discriminates 
form instance a cc pair with the binding energy of a T' but the radius of a Xc and so 
forth. For the J/T the minimum radius tq = 0.1 fm eliminates too high values of the 
binding energy. Such requirements do not apply to the T , the binding energy being in this 
case limited by the depth of the potential (controlled by the cutoff A, as discussed in the 
previous subsection). 

The time evolutions of the average size of pairs thus prepared are presented in Fig. 5 
for different temperatures. The harmonic oscillator pattern expected from the analysis of 
Sect. 4.3 for pairs of small initial sizes is clearly visible. There are indeed cases where 
(Tqq) remains almost constant for a certain time interval, reflecting the fact that the cor¬ 
responding pair is highly correlated, or “bound”. The lower the temperature, the longer the 
correlation lasts. One also observes the expected “sequential melting” of T', Xc J/^ and 
T as temperature grows. Of course the sequential dissociation of T', Xc J/^ just reflects 
the inequalities of their respective sizes, r^i > > TJj^^. One may try and attribute dif¬ 

ferent “melting temperatures” to the dissociation of the various bound states. For example, 
we see from Fig. 5 that the initial plateau associated with the average T' radius is absent 
at T = 220 MeV, indicating that immediately dissolves at this temperature, while the 
plateau is still visible at 190 MeV. One may then infer that the melting temperature of T' is 


^^The binding energy is known at each time step of the simulation, since we follow both the velocities 
and the positions of the particles. 
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Figure 5: (Color online.) Average quarkonia radius as a function of time. The pairs are prepared as “bound 
states” following the procedure explained in the text. The shaded part indicates the region in which the 
charmonia radii are smaller than the Debye radius vd- Statistical errors are too small to be plotted. In 
each plot, the upper curve represent the less bound system, the lower curve the most bound one. 


T ~ 200 MeV. Using the same argument of the size of the screening radius, we can extract a 
melting temperature of T ~ 310 MeV for Xc ■ On the other hand, it is evident that the J/^ 
survives up to much higher temperatures than the other two charmonium states. However, 
for the J/'h we can not use the above strategy to estimate its melting temperature, because 
of the limitation of the numerical setup: when the temperature increases {T > 400 MeV) it 
becomes impossible (with the present choice of parameters) to prepare an initial J/^ with 
AE > 550 MeV, the potential well is simply not deep enough (see Fig. 2). Later, we shall 
estimate the melting temperature of the J/'ll by using a different procedure. 

In the last panel of Fig. 5 we also compare the Xc, J/^ and T behaviours at T = 280 
MeV. We see that the average bb pair is far more strongly correlated than the cc pair, and 
the T radius remains small ((rx ) < J’d ) for a relatively long time (we shall see in the next 
subsection that the melting temperature of the T (IS) state is T > 600 MeV). 
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T (MeV) 

160 

190 

220 

280 

Xc 

(40-43)% 

(28-30)% 

(16-17)% 

(1-2)% 

T' 

(12-14)% 

(7-8)% 

0 % 

/ 


Table 1: Fractions of Xc and 4*' eventually becoming J/^’s. 


The curves representing the Xc exhibit an interesting phenomenon. One sees that in 
all cases, the corresponding radius tends to decrease initially, bringing the Xc closer to a 
more stable bound state (J/T). This is a clear indication that, at these temperatures, the 
binding forces are not yet entirely screened. While on average, the relative kinetic energy 
prevents the Xc to really decay into a J/'k, as the curves in Fig. 5 indicate, a substantial 
fraction of the pairs prepared as Xc does decay into J/ip's, as shown in Fig. 6. It is possible 
to estimate the percentages of Xc and 'll' states that decay into J/T, a process known as 
feed-down. In Table 1 the feed-down percentages of Xc and T' are listed for some values 
of the temperature. We found that, for each temperature, there are more Xc than T' 
states that decay into J/'k , and the feed-down mechanism decreases when the temperature 
grows: the more fragile states prefer to dissociate rather than form a more strongly bound 
system. Amusingly, the feed-down fractions obtained here at T = 190 MeV are similar to 
the experimental values quoted in [50], although of course the physical context is rather 
different. In Fig. 6 we compare the different behaviors of the Xc and T' average radii, 
separating those which decay from those which do not. Looking on the left of Fig. 6 we 
notice that even the non-decaying Xc states initially reduce (on average) their radius (also 
at T = 280 MeV, as seen in the last panel of Fig. 5). This why their average lifetime (see 
Table 2) remains almost the same below the melting temperature (~ 310 MeV), whereas 
the average lifetime of a non-decaying diminishes as the temperature goes up. 

The lifetimes of Table 2 are calculated by averaging the time intervals needed for the 
radii of Xc 'k' (both not-decaying) and J/'k to become larger than the Debye screening 
length. One notices that the J/'k lifetime at T = 280 MeV is still quite appreciable. 
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Figure 6: (Color online.) On the left: Comparison between (rcc) of Xc states that become J/^ (« 41% - 
lower curve) and those that do not decay (« 59% - upper curve). On the right: Same comparison for . 
The pairs are prepared as bound states as indicated in the text. For similar initial conditions, they evolve 
statistically to different final states. 


T (MeV) 

160 

190 

220 

280 


>10 

>10 

4.9± 0.2 

2 .8± 0.2 

Xc 

1 .6± 0.1 

1 .6± 0.1 

1.5± 0.1 

1 .6± 0.1 


0.7± 0.1 

0.5± 0.1 

0 .1± 0.1 

0 


Table 2: Average charmonium lifetimes (in fm/c) in the quark-gluon plasma. Only the Xc and that do 
not decay into J/T are taken into account in the lifetime estimates. 
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5.3. Many heavy quark-antiquark pairs 


In the previous subsection, we saw how a single heavy quark antiquark pair can evolve 
from an apparent bound state to a system of two independent quarks that eventually ther- 
malize with the plasma on long time scales. We could also observe, with a proper selection 
of the initial conditions the expected phenomenon of sequential dissociations. Finally, we 
provided some criterion to get a crude estimate of the lifetime of the bound state in the 
plasma. We would like now to examine how these features are modified when several 
pairs are present in the plasma. The simulations that we shall present were performed for 



_^_I_^^_I_1_ 

160 180 200 220 240 260 280 

T (MeV) 


Figure 7: Average total energy of a system of 10 cc pairs in thermal equilibrium as a function of temperature. 
Before measuring the energies, we ran the simulations for a time interval of 100 fm in order to let the system 
thermalize. Simulations were performed in a periodic cubic box of side 4 fm and statistical errors are again 
too small to be plotted. 


= 2,10 and 50 quark-antiquark pairs, in a cubic box of side 4 fm, with periodic boundary 
conditions. 

When there are enough pairs in the system, one expects them to evolve towards an ideal 
gas of the constituents, if the temperature is high enough. The average energy for a system 
of = 10 pairs is plotted in Fig. 7, and compared to that of an ideal monoatomic gas of 
2N particles, 

E^^ = ^{2N)KbT. (5.81) 
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The expected trend is clearly visible, and at the largest temperatures considered, T > 280 
MeV, the ideal gas limit is almost reached. At such high temperatures, most of the pairs 
dissociate if one waits long enough. On the other hand, at lower temperatures, pairs may 
survive and this results in the average energy of the system being lower than that of the ideal 
gas at the same temperature. We note that the process of dissociation, considered from 
this thermodynamical point of view, is a gradual process: even at high temperature there 
remains some finite probability to find a bound pair. Given the length of the simulation 
(over 100 fm/c), and that, in this range of temperatures, a single pair would eventually 
dissociate, the equilibrium state that we are observing results from the balance of the two 
competing effects of dissociation and recombination, as we shall discuss in more details 
shortly. 



0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 


distance (fm) 


Figure 8: The figure shows the distribution of distances to the nearest antiquark from a given quark. This 
probability is computed in the following way: once thermal equilibrium is reached, one computes from 
each quark the distance to the nearest antiquark, draw an histogram, and normalize in order to get the 
distribution. Simulations were performed for a system of 10 cc pairs in a cubic box of side 4 fm, with 
periodic boundary conditions. 


The presence of bound pairs in the system can also be inferred form the analysis of 
another quantity that is directly sensitive to the correlations between two particles, namely 
the probability distribution Pqq of the distance from a given quark to the nearest antiquark. 
In an ideal gas, this distribution is given by 


pideal 

^qq 



iV>l 


3 

a 
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where a 



is the mean distance between the antiquarks and p 


y the density 


of antiquarks. The peak of the ideal gas distribution for = 10 quark-antiquark pairs in a 
cubic box of side 4 fm, is at rpeak = ~ 1-15 fm. This peak is clearly visible in the 

distribution Pqq of the interacting system which is plotted in Fig. 8. But this figure reveals 
also another feature: at low temperature, there is also a sharper peak reflecting the presence 
of highly correlated states in the system. These, we associate with the bound states. In line 
with the previous plot. Fig. 7, this peak disappears when T > 280 MeV. From Fig. 8 we can 


also infer that a correlated cc pair has a maximum radius of approximately 0.3 fm, which 
is indeed similar to the values of the Debye screening length in this range of temperature. 
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Figure 9: On the left: Average recombination time as a function of temperature. On the right: Probability 
of recombination times for three values of temperatures. The quark-antiquark distance has been chosen to 
be less or equal to 0.3 fm for the quark-antiquark configuration to be considered a pair (see text). Both 
simulations were performed with a system of 10 cc pairs in a cubic box of side 4 fm, with periodic boundary 
conditions. 


We turn now to a more detailed study of the process of recombination. We start with 
the evaluation of the recombination times for the pairs as a function of the temperature, 
that is the average time needed for a quark (antiquark) to form a pair, once the quark 
(antiquark) moves away from its previous antiquark (quark) partner. In doing this cal¬ 
culation we carefully avoid counting the contributions of “non-interacting” events, that is 
the occurrences where a quark passes by an antiquark without forming an actual pair. In 
order to eliminate such events, we performed simulations for a non-interacting system with 
a constant (space-independent) drag constant (see Eq. (4.71)) and we calculated the corre¬ 
sponding normalized distribution Pfree(t) of the time intervals t in which a charm and an 
anticharm stay close together within a sphere of radius 0.3 fm. Then, for each temperature, 
we define a minimum lifetime r by the condition 

dt Pireeit) < 1% . (5.82) 
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By selecting pairs that stay together for a time greater than r, only pairs formed because 
of the interactions (and not those resulting from random encounters) contribute to the 
recombination times. Note that the procedure does not allow for a detailed analysis in 
terms of various bound states, as we were able to do for the dissociation: the small lifetimes 
typical of Xc and 'h' are automatically discarded by the procedure, so that we implicitly 
consider only J/^ regeneration. 

As one can see on the left of Fig. 9, the outcome for a system of 10 cc pairs is that the 
recombination time increases almost linearly with the temperature, starting from a value of 
tree = (62.9±2.5) fm at T = 160 MeV and reaching a value of tree = (185.5 ±6.8) fm at T = 
280 MeV. As one increases the temperature one increases the rate of encounters, but also the 
relative kinetic energies of the pairs, preventing binding. Another important observation 
is that the recombination times are very long, so long that one may wonder whether the 
mechanism of recombination could be of any phenomenological relevance. However, as the 
graph on the right panel of Fig. 9 shows, the distribution of the recombination times is 
very broad. Thus, even if the standard errors of the graph on the left panel of Fig. 9 are 
small (because of the large statistics), the standard deviations are of the same order as the 
average values: over the lifetime of the quark-gluon plasma (~ 10 fm/c) there is effectively 
no characteristic time scale for recombination. 

One can nevertheless push the discussion a bit further and quantify the process in a 
simple way. Note first that the recombination time is expected to go up when the number 
of particles decreases. This is indeed what we obtain from our simulations. We find that 
the average recombination time is, to a good approximation, inversely proportional to the 
number of pairs present in the system: trecA^ ~ , with A/j a (temperature-dependent) 

recombination rate. This effect is also (qualitatively) visible in Fig. 10 that displays the 
fraction of surviving J/T (and T) particles as a function of time, for different number of 
pairs in the system: one notices that recombination events are more frequent in a system 
with a greater number of cc pairs. One may also observe that the effect of recombination 
becomes relatively more important as the temperature grows. This is visible for instance 
from the development of a plateau suggestive of equilibrium that is most clearly seen at 
the highest temperature (T = 220 MeV). Finally, the last panel of Fig. 10 compares the 
behaviors of cc and bb at a given temperature over a long time scale. One sees that there is 
a lapse of time before the bb bound state starts to “feel” the action of the thermal medium. 
This time delay to is about to ~ 4 fm/c. A similar effect also occurs for charmonium, but for 
a smaller to < 1 fm/c. This dependence on the mass is a clear indication of the important 
role of the collisions in the dissociation process. 

At the same time, the effect of the binding forces is certainly also present. This we see 
indirectly by studying the cutoff dependence of the results. To that aim, we have repeated 
simulations for various values of the cutoff. As we have seen earlier, the dominant effect 
of a change in the cutoff is to change the depth of the potential. A larger cutoff leads to 
a deeper potential, and a longer lifetime, and this effect persists up to values of the order 
A ~ 6 GeV, above which it attenuates considerably. In turns, this alters the recombination 
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rates since the pairs with too short lifetimes are eliminated by the procedure with which 
we identify bound pairs. 

One may understand quantitatively the behaviors identified in Fig. 10 from a simple rate 
equation. Let us denote by X£){T) the dissociation rate and by Xji{T) the recombination 
rate. Both are functions of the temperature. The rate equation describing the time evolution 


T = 160 MeV 


T= 190 MeV 



2 initial pairs 
10 initial pairs 
50 initial pairs 


time (fm/c) 


time (fm/c) 


T = 220 MeV 


T = 190 MeV, 10 initial pairs 




Figure 10: Fraction of surviving pairs as a function of time for three different temperatures. The N pairs 
{N = 2,10, 50) are prepared at t = 0 as bound states as indicate earlier in the text. The short time behavior 
is dominated by dissociation. The process starts however only after some small delay to ^ 1 fm/c. This 
delay is much longer the the bottomonium, as revealed by the comparison displayed in the bottom-right 
panel: fraction of surviving J/'F and T, for a system of = 10 c-c or b-b pairs. Simulations were performed 
in a periodic cubic box of side 4 fm. 


of the number of surviving QQ pairs N{t) is (see also [10]) 


= -XoNit) + XRN,{t)Ng{t) , 


(5.83) 


where Nq = Ng is the number of free heavy quarks (or antiquarks) in the plasma. Equation 
(5.83) together with the initial condition N{t = to) = Nq , Nq(tQ) = A^q(to) = 0 can be 
analytically solved for the fraction of surviving pairs: 


N{t) _ 1- ^tanh(f(f-fo)) 

Nq 1 + ^ tanh (^(f - fo)) ’ ~ ° 


(5.84) 
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Figure 11: Left panel: fit of the solution (5.84) of the rate equation for a system of 10 initial J/'i> at 
three different temperatures. We notice that eq.(5.84) fits the curves very well with the recombination 
rates Ah equal to the ones derived in Fig. 9. Right panel: high temperature behaviour of charmonium and 
bottomonium. 


where = \/Xd{^d + 4A/jA^o) and we used Nq{t) = Nq{t) = Nq — N{t). The time to is the 
time at which dissociation starts, as defined earlier, and this time is chosen as the initial 
time when solving the rate equation (5.83). 


T (MeV) 

V (Wc) 

(fm/c) 

160 

23 

625 

190 

9.2 

1000 

220 

4.6 

1350 


Table 3: The inverse of the dissociation rate A^^ and the recombination rate Ajj^, for various temperatures. 


From the results of the simulations, and using Eq. (5.84), we can extract values for 
the dissociation and regeneration rates. The results are reported in Table 5.3 for a few 
temperatures. The quality of the fit can be seen on Fig. 11. The values of the dissociation 
rate thus determined can be used to infer the lifetimes of the bound states. The values 
of the J/T lifetimes obtained from the present fit agree with those previously obtained by 
analyzing the time evolution of the J/T radius (see Table 2). The analysis of the values 
of the dissociation rate just obtained suggests that at a temperature T = 600 MeV, the 
lifetime of the J/il) is still non vanishing, and is about .5 fm/c. At the same temperature, 
the lifetime of the T is about 1.5 fm/c. The fit of bottomonium data gives essentially zero 
recombination rates already at T = 160 MeV (Ajj ~ (1 ± 1) • 10“® fm/c) and much smaller 
values of the dissociation rate as compared to the charmonium case. These numbers reflect 
of course the greater stability of the T as compared to the J/T. 
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6. Conclusions 


We have presented an approach that treats in a unique framework most of the important 
aspects of the evolution of a collection of heavy quark-antiquark pairs propagating through 
a quark gluon plasma. The approach starts from first principles, and leads, through well 
defined approximations, to a complete dynamical description with a unified perspective on 
many different physical effects, usually treated with different models. Of course, several 
approximations are needed to arrive at tractable calculations. However these approxima¬ 
tions can be improved, and their presence should not obscure the overall consistency of the 
general scheme. We find it particularly important, for instance, in view of their relevance 
for the interpretation of the data, to have the processes of dissociation and recombination 
treated on the same footing. 

The main question that is addressed in this paper is of a general nature, it concerns 
the fate of a collection of heavy quark pairs in a hot and dense environment, with the 
possibility for these heavy particles to form bound states. This does not involve QCD 
dynamics in an essential way, and this is the main reason why we have restricted ourselves 
in this paper to Abelian plasmas where the same question can be addressed in a much 
simpler setting. Specific features of QCD can be implemented within the present scheme, 
with perhaps some approximations becoming less accurate. In particular, one of the main 
approximations can be understood as a weak coupling approximation, which consists in 
neglecting the non linear coupling of the gauge (Coulomb) field with which the particles 
interact. The field fluctuations are then Gaussian, which allows for a simple calculation of 
the influence functional in terms of a 2-point function that characterizes entirely the effect 
of the plasma on the heavy particles. In QCD, the non linear couplings are not as strongly 
suppressed as in QED, and the approximation may be less accurate. 

Further approximations lead to a classical treatment of the dynamics in terms of a 
Langevin equation, in which the noise term accounts for the effect of the collisions between 
the heavy particles and the plasma constituents. This noise terms depends on the positions 
of the heavy quarks at each time steps. This dependence is an important aspect of the 
dynamics. 

The simulations presented in this paper are the first in a program that can be improved 
in many ways. Some of the approximations that have been made can easily be relaxed, 
such as for instance the Abelian approximation, as we have discussed already. The classical 
treatment of the dynamics through a Langevin equation could be improved, e.g. including 
the leading order quantum corrections to the Langevin equation, as shown in [57]. The 
basic ingredients such as the transport coefficients, can be calculated with greater accuracy. 
Finally, once some of these improvements are implemented, more realistic phenomenological 
applications can be envisaged. We hope to be able to report on some of these developments 
soon. 
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Appendix A. Taylor expansion of the Influence Functional 

In this Appendix, we perform the Taylor expansion of $[<5] = ^qq[Q] + + 

obtained from Eqs. (4.57) and (4.58) of section 4 after performing the change of 
coordinates of Eq. (4.59). We first analyze the contribution 4 *qq[Q] involving only the 
heavy quarks. We have 



(A.l) 


where the time dependence is hidden in the coordinates R and Y. The dot symbol in this 
expression, as well as in the rest of this section, denotes a scalar product and involves the 
three cartesian components of the vectors. We want to expand the expression (A.l) to 
second order \n y. To do so, we use the well-known Taylor expansion of a scalar function 
/ of a n-dimensional vector 

f{x) = /(a) -F (® - a) • V/(a) + ]^ {x - a) ■ {a) ■ {x - a) + ... 

where x = (xi,..., Xn), and the gradient and Hessian matrix are given, as usual, by 


V«/(a) 




d^fjx) 

dxadxp 


aj=a 


a, /3 = 1, • • • , n. (A.2) 


By applying this formula to the hrst line of Eq. (A.l) we get (with Yji = rj—r^, y ^ = Vj—yj) 


^(rj* - \yji) -V{rji + Ivji) = -Vji ■ VK(rjj). 




(A.3) 
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Similarly 

-W(rji - - W(rji + = -2W(rji) - ^ • n{rji) • y^^, (A.4) 

and (with = pj + yj) 

W (r,, + ^) = W{r,,) + ^ ^ ^ (A.5) 


Note that the middle term in the right hand side of the last equation will disappear in the 
summation over i and j, since it is antisymmetric {'VW{rji) = — V14^(rjj))^^. 

Let us now consider the terms that involve the time derivative. We write this as 


2 


■ji+ 2 


d 

drji 


W 


■j* 


+ nV 


J* 


(A.6) 


and use the expansion of W above. When keeping only the symmetric terms, i.e., those 
which survive in the summation over i and j, this yields 



• ^ • Vji + i/ji • VW . 


(A.7) 


At this point, we note that one can write pj^ • VW (rji) as —pji ■%{rji) • rji after integrating 
by part in the integral over time appearing in ^qq ■ The boundary terms coming from this 
integration by parts vanish because the coordinates Qi and Q 2 coincide at both ends of 
the Schwinger-Keldysh contour, that is 


= QjAtf) - fljAtf) = 0, yj{ti) = Qj,i(ti) - qjAti) = 0- 


Collecting all intermediate results, we get 

2 N 

foolR.YI = / tid-ra.i- 

+ (3 (pji-nirij) -rji- 


VC {vij) + i {pji • n • Pji - Pji 
Pji-nirji) ■Vji'^ , 


■ ^ ■ yji) + 

(A.8) 


which we can rewrite as 

fl2 ^ rtf 

= -^Y^ / dt[2yi-Wirij)--iyi-nirij)-yj + ^yi-'H{rij)-rj] . 

(A.9) 


^^The Hessian matrix of W is the only such matrix in the present discussion, so we denote it simply by 
H, without any explicit reference to W in the notation. That is, in the notation of Eq. (A.2), H = . 

^^When using regularized potentials L(r) ad W{r) such that VL(r) and VlT(r) both vanish at r = 0, 
the same cancellation holds for the terms with i = j. 
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The result for is obtained trivially from ^qq via the change of variables y ^ y,r ^ r. 
Let us then consider the expansion of the remaining term, Y], We have 



(A.IO) 


with now Tji = Tj — Fj, = y^ — y^ and y^^ = y^ + y^. By using similar manipulations as 
above, one finds that the last two terms contribute 

^ {Vi • VW(rjj) - ii • n{rji) • y^-^) = (y^ • n{rji) • Vji + f* • n{vji) • y^-^) , 

(ATI) 

where we have used an integration by part. Moving up to the second line, we get 

i W (^rji + ^yji)^ + i W (^rji - ^yji)}^ = 2iW (rji) + ^yj^ • n(rji) • yj^ (A.12) 

and 

-i W (^rji + ^yji)^ - i W (^rji - = -2i W(rji) - ^yj^ ■ n(rji) ■ y^^ (A.13) 

As for the first line, it yields simply 

V (^rji - ^yji)^ - V (^rji + ^yji)^ = -yji ■ W(rji). (A.14) 

Collecting all the intermediate results, we can then rewrite the influence functional as follows 


^ rtf 

^pp[R->Y] = dt{-{yj-yf)-VVirj-fi)+ iyj-nirj-Ti)-yf 

i,j=i "'L 

- rj) ■ Tj + yj ■ nirj - r*) -fi)} . 

(A.15) 


By collecting the Taylor expansions of ^qq , ^qq derived in this appendix, and 

using the definitions (4.63) and (4.64), one easily obtains the equations (4.61, 4.62) of the 
main text 
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Appendix B. Derivation of the generalized Langevin equation 

In this Appendix, we show that the dynamics encoded in the path integral (4.61) is 
equivalent to that described by the generalized Langevin equation (4.66). Let us start by 
considering the Langevin equation for a particle of mass M moving in an A-dimensional 
space: 


M n =-Mjijf-j + fi{r) + (i{t) , i = l,...,iV, (B.l) 

where r* denotes a coordinate of the particle, fi and fj its first and second time derivatives, 
fi is an external deterministic force, and a white stochastic force with the following 
properties 

{m)i = 0, (C*(t) 6(0 )? = Ay5(t-0, A., =2MT7,,-. (B.2) 

Here 7 is a real symmetric matrix, and we have used Einstein’s relation between the noise 
and the dissipative terms. The equation that we need to consider is a generalization of 
Eq. (B.l) in which the matrix 7*^ (and hence Xij) depends on the position r of the heavy 
particle. It is of the form 


M fi = -M'yij{r) fj + /^(r) + 6(r, t ). 


with a so-called multiplicative noise 


(B.3) 


6(r,t) :=Wij{r) 6(6, Wik{r)wjk{r) = Xij{r), { Ci(t) ij{t') )^ = 5ij S{t - t'). 

(B.4) 


An equation such as Eq. (B.3) may suffer from discretization ambiguities in the case where 
the inertia term, the left hand side of the equation, is ignored, leading to the so-called 
“overdamped” Langevin equation (see for instance [51]). These ambiguities reside in the 
choice of the point r where the noise is evaluated when one solves the stochastic equation. 
One may indeed choose to evaluate the noise w{r) at any point r between r{t) and r{t + At), 
where At is the discrete time step. This leads to an uncertainty of order 


drc 

dr 


• f At 6 


(B.5) 


with w{r) assumed to be a smooth function of r. In the overdamped case, we have r ~ ^ ~ 
l/\/At so that the uncertainty is of order unity and remains hnite as At —)■ 0. However, 
as discussed in [52], such ambiguities may not appear when the inertial term is present, 
which is the case of interest in the present discussion. This is because, one can rewrite the 
equation (B.3) as a set of two coupled equations. 


f = t; 

Mi} = —M'y{r) ■ v + f(r) -|- m(r) • ^(t). (B.6) 
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In this case, while h ~ ~ v itself remains hnite, and so does f. It follows that the 

uncertainty (B.5) is now of order '/At and it vanishes as At —)■ 0. 

We shall then proceed to the discretization of the system of equations (B.6), and to be 
specihc, we shall use the Ito convention, where the noise is estimated at the position of the 
particle before the time step considered. The discretized form of Eqs. (B.6) reads then 

r^n 

rW _ r(^-i) = At + At / ds C{s) 

1 

M + Atf^"”^) + T” ds , 

(B.7) 

with initial conditions = tq and = Vq. In the equations above, we have set 
^(n-i) _ and in order to simplify the 

notation. We have 

(B.8) 

To write Eqs. (B.7), we have divided the time interval [0,t] into n time step At, At ;= 
tn — tn-i , n = 1,..., h. We have also introduced in the hrst equation an auxiliary white 
Gaussian noise Q{t) with properties (see also [53] for a similar procedure) 

( Ci{t) ) = 0, ( C,i{t) Cj/) ) = ti 5ij 5{t - t'). (B.9) 

The additional factor At in front of the integral of the noise in the hrst equation (B.7) 
ensures that r = t; in the limit At —)• 0, thus avoiding any discretization ambiguity. 

The conditional probability for the particle to be found in the conhguration (r^""), 
at time tn, given that it is in the conhguration at time tn-i, can be written 

as [51] 

P (rW,^W,tn I r(-i),^(-i),tn-i) = (5 (rW - (<5 )«> 

(B.IO) 

where = rsoi(,tn;/'^~^\tn-i) and = Vsoi{tn;v^^~^\tn-i) are the solutions of the 
discretized equations (B.7) for a given realization of the (independent) noises ^ and (), and 
given values of and at time tn-i. The probability of having a hnal conhguration 

^■^In [53] the extra At is not used, but fi is eventually sent to zero. The advantage of keeping the factor 
At explicit is that it makes obvious that all potential discretization ambiguities disappear as At —>■ 0. 
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(r("'),v(”)) = (r, v) at time tn = t given an initial configuration = (ro,t7o) at 

time to = 0 is given in terms of the probability (B.IO) of an elementary step by 


P(r, v,t I ro, vo,to) = 


n—1 „ „ 

n=l 


This formula is the starting point for building the path integral. In order to evaluate the 
average over the noises of the delta functions in Eq. (B.IO), it is convenient to define the 
following A^-dimensional vectors 

pin 

g{n) ^(n) _ ^(n-1) _ ^(n-1) _ (^(t) (B.12) 

'^^n —1 

pin 

;= M + MAt - At - dt 

tn — 1 


These functions vanish respectively when and and we have 








(B.13) 


In order to calculate the noise averages, we use the Fourier representation of these delta 
functions. We get then Taking the mean value of these quantities, we get 





2tx 


e * 








' —OQ 

f+co 




27r 


(B.14) 


where we have exploited the presence of the explicit factor At in order to evaluate the 
average of the last exponential factor to linear order in At, where it reduces to unity. 
Similarly, 



X 


f + OO 


(n) 




(e 


. (n 

'Vk 


dt/fc 

27r 






(B.15) 


where now we need to push the expansion of the last exponential to second order 




ds 


' ^n —1 


' tn — ^ 




1 1 (n) (n) Pn-l) . , -i At 

= 1 - Vi Va ■ At PS e 2 yfc yj kj 


2 yj ''kj 


(B.16) 
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The probability P(r, v,t | rg, vo,to) of EQ- (B-11) can therefore be written as (to within 
the factor coming from the delta functions (B.13) and that can be absorbed in the 
normalization) 


n—1 






n=l ' 


X exp 


• —oo J —oo 



1 

1 

1 

1 

_1 


V 


( \ I — V 

■ M 




(B.17) 


At this point we note that we can integrate over 2 , thereby reconstructing the delta function 
(5 [r(t) — which removes the integration over the velocity v. Note that there is no 

integration over Zq nor Zn, in line with the fact that Vq and Vn are fixed by the initial 
and final conditions on the paths. Note also that the relevant path are differentiable, with 
r = V finite. After integrating over Vq and v, and sending At —)■ 0, the remaining path 
integrals over r and y yield the probability P{rf,tf \ Tipi) in the form 


P{rf,tf I Tipi) = j Dr j Dy exp 


rtf 


\_J ti 


dt £(r,y) 


(B.18) 


with 


>C(r,2/) 


—iy • (Mr + M^{r) • f — f(r)) 


1 

2 


y • A(r) • y. 


(B.19) 


where we should remember that there is no integration on the end point of the y path 
integral, or equivalently that y(t/) = y(tj) = 0. The structure of this expression is identical 
to that of the conditional probability (4.61) derived in section 4.2. 


Appendix C. Numerical algorithm 

In order to perform numerical simulations we use an explicit second-order algorithm^^ 
which requires the evaluation of a single function at each time step. This algorithm can be 
summarized as follows. At each time step: 

• use an orthogonal transformation to pass to coordinates for which the real and sym¬ 
metric matrix 7 is diagonal; 

• perform the stochastic Verlet algorithm (see below); 


^®By second-order algorithm we mean that the convergence of the algorithm is of the order of Ap, with 
At the time step of the simulation. 
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• come back to original coordinates. 

Let us detail the second step. After diagonalizing 7 , we are left with with N independent 
stochastic equations of the Ornstein-Uhlenbeck type [54] 

P = -li'r) P + f{r) + y^2MTj{r) ^{t) , (C.l) 


with 



(C.2) 


{m) = o, = (c.3) 

Eq. (C.l) can be written as 

dp = — 7 (r)p dt + /(r) dt + ^/2MT 7 (r) (C.4) 

which admits the exact solution [54] 

Pt+At=Pie-“ + ^(l-e-“) + yiWT0:^e=^6, (C.5) 


where a(r) = 'y{r)At. This result enables us to write the propagation of Eqs. (C.l) and 
(C.2) in a manner similar to the Verlet algorithm used in Newtonian dynamics [55] 


r+ 


Pt+At 

^i+At 


n + 


Pt 

2M 


At 


(1 - a(r +))Pt + /(r+)At + ^/2MTaJj^it 
Pt+At 


r+ + 


2M 


At. 


which is often referred to as the stochastic Verlet algorithm [56]. 


(C. 6 ) 
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